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Abstract 

We study effective shear viscosity /i* and effective extensional viscosity A"^ of concentrated 
non-colloidal suspensions of rigid spherical particles. The focus is on the spatially disordered 
arrays, but periodic arrays are considered as well. We use recently developed discrete network 
approximation techniques to obtain asymptotic formulas for ^* and A* as the typical inter- 
particle distance 5 tends to zero, assuming that the fluid flow is governed by Stokes equations. 
For disordered arrays, the volume fraction alone does not determine the effective viscosity. Use 
of the network approximation allows us to study the dependence of /i* and A* on variable 
distances between neighboring particles in such arrays. 

Our analysis, carried out for a two-dimensional model, can be characterized as global because 
it goes beyond the local analysis of flow between two particles and takes into account hydro- 
dynamical interactions in the entire particle array. Previously, asymptotic formulas for ^* and 
A* were obtained via asymptotic analysis of lubrication effects in a single thin gap between two 
closely spaced particles. The principal conclusion in the paper is that, in general, asymptotic 
formulas for /i* and A* obtained by global analysis are different from the formulas obtained 
from local analysis. In particular, we show that the leading term in the asymptotics of /i* is of 
lower order than suggested by the local analysis (weak blow up) , while the order of the leading 
term in the asymptotics of A* depends on the geometry of the particle array (either weak or 
strong blow up). We obtain geometric conditions on a random particle array under which the 
asymptotic order of A* coincides with the order of the local dissipation in a gap between two 
neighboring particles, and show that these conditions are generic. We also provide an example 
of a uniformly closely packed particle array for which the leading term in the asymptotics of A* 
degenerates (weak blow up). 
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1 Introduction 



Concentrated suspensions are important in many industrial applications such as drilling, water- 
coal slurries transport, food processing, cosmetics and ceramics manufacture. In nature, flows of 
concentrated suspensions appear as mud slides, lava flows and soils liquefied by the earthquake- 
induced vibrations (0, [Z]) [TB]). 

An asymptotic formula for the effective viscosity of a suspension of non-colloidal particles in a 
Newtonian fluid, derived in [S], is based on the local lubrication analysis of the energy dissipation 
rate in the narrow gap between a pair of nearly touching particles. The distance between two 
neighboring particles in a periodic array is the small parameter in the problem. For periodic 
arrays, this inter-particle distance is uniquely determined by the volume fraction of particles, so 
that the asymptotics of the effective viscosity is obtained as a function of the volume fraction (p 
that is close to the maximal packing volume fraction (j)rcp- the asymptotics of the effective viscosity 
obtained in [H] has the form 

Ae-^ + 0{lne), (1.1) 

as e ^ 0, where e = 1 — {4>/ 4>rcpY^'^ ■ The formulas for effective viscosity of periodic suspensions 
in the whole space (without boundary) subject to a prescribed linear flow, obtained in [TT] . 
also rely on the local lubrication analysis. Asymptotic representations for the components of the 
effective viscosity tensor calculated by |11| are of the form 

Ae-^ + B\ne + 0{1), (1.2) 

Recently, concentrated random suspensions were investigated numerically by Jl] using accelerated 
Stokesian dynamics. It was observed that the behaviour of the effective high frequency dynamic 
shear viscosity of disordered suspensions can be accurately described by the asymptotic Sine, 
indicating degeneration of the leading term in the asymptotic expansions 1)1.2(1 (weak blow up). 
The authors of also show that their results are in good agreement with available experimental 
data (Cni) \n\)- This suggests that for generic random suspensions, the asymptotics of the effective 
viscosity defined by the (properly normalized) global dissipation rate cannot be identified with the 
local dissipation rate in a single gap. 

In this paper, we use the discrete network approximation proposed in |4j to study the asymp- 
totics of the shear effective viscosity and the extensional effective viscosity A* corresponding 
to general disordered particle arrays. For such arrays, the volume fraction alone is not sufficient 
for determining the effective viscosity. Therefore, instead of e, we use the inter-particle distance 
parameter 6 that controls the distances 5ij between neighboring particles. 

The mathematical construction of [3] accounts for the long range hydrodynamical interactions 
between the particles, and provides an algorithm for calculation of the effective viscosity, which takes 
into account variable distances between neighboring particles in non-periodic arrays. Furthermore, 
in ^ it was observed that the leading term of the asymptotics may degenerate due to the external 
boundary conditions and geometry of the particle array, while in the scalar case (|S]) the order of 
the leading term is the same for all particle arrays that form a connected network. This paper is 
devoted to a detailed study of this degeneration phenomenon. In particular, we clarify the issue of 
weak versus strong blow up in the asymptotics of the effective viscosity. 

The definition of the effective viscosity employed in this paper differs from the definition in |llj . 
Their definition is closely related to the definition of the effective material properties of periodic 
media in homogenization theory (see e.g. 0, jlUj). where these properties are defined via analysis 
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on a periodicity cell. Then these properties are determined by the properties of constituents and 
geometry of the periodic array, independently of the external boundary conditions and applied 
forces. Such properties are usually defined in the limit when particle size tends to zero while the 
number of particles approaches infinity. 

Our definition is more directly linked to the viscomctric measurements that necessarily involve 
boundary conditions. We assume that particles of a fixed size are placed in a bounded domain 
(apparatus), and the typical distance between the neighboring particles approaches zero. Thus, 
the total volume fraction of particles approaches the maximal packing fraction. In this case the 
effective viscosity will be influenced not just by the inter-particle interactions but also by the forces 
or velocities prescribed at the boundary of the apparatus. Since we work with linear models, 

(A*) are independent of the applied shear (extension) rate. However, we also show that the 
asymptotic order of the effective viscosities depends on the relation between the orientation of the 
velocity prescribed at the boundary, shape of the boundary, and orientations of the line segments 
connecting pairs of neighboring particles. 

We present sufficient conditions for degeneration (weak blow up) and non-degeneration (strong 
blow up) of the leading term. We show that /x* always exhibits weak blow up, while local analysis 
alone predicts strong blow up. We also show that the asymptotic order of A* depends on the 
geometry of the particle array. In the paper, we define a broad class of arrays of particles for 
which the leading order term of A* does not degenerate. This situation is typical in the sense that 
for a generic array /x* and A* have vastly different values, and their ratio depends on the inter- 
particle distance, which indicates possible non-Newtonian behaviour of the effective fluid under the 
imposed boundary conditions. Wc also give an example of a closely packed particle array for which 
the leading term in the extensional viscosity degenerates. 

2 Mathematical model 

We consider a concentrated suspension of rigid, non-Brownian, neutrally buoyant particles in a 
viscous incompressible Newtonian fluid. In the two-dimensional model, the suspension occupies 
a domain CI which is a square of side length two centered at the origin. The boundary of Q is 
denoted by dfl. The upper and lower sides of dfl are denoted by dU'^ = {x. : X2 = 1} and 
di}~ = {x : X2 = —1}, respectively. We also let ei,e2 denote the Cartesian basis vectors parallel 
to the sides of dQ. The particles , j = 1,2,...,N are modelled as rigid disks with centers x-?, 
placed in Q,. For simplicity, we consider the monodisperse case, so that all particles have radius a. 
The part of CI which is not occupied by particles is the fluid domain, denoted by Ctp- 
The fluid at low Reynolds number is governed by Stokes equations 



where iJ, is the fluid viscosity, v the velocity field, and P is the pressure. 

In this paper, we consider the external boundary conditions of the shear and extensional types. 
The shear type boundary conditions are given by 



/xAv-VP = 0, divv = 0, inftF- 



(2.1) 




on dQ~^, 
on dO,~ , 



(2.2) 



where 7 is a constant shear rate. 
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In the case of extensional boundary conditions the velocity is prescribed as 

e{—e2 + xiei) on dfi^, 



(2.3) 

e(e2 + xiei) on (90 , 

where e is a constant extension rate. The two remaining (vertical) sides of 9$7 are free surfaces, 
where the zero traction condition is prescribed. 

To define particle velocities, we first recall that a rigid body moving in the plane defined by 
ei,e2 has a velocity vector of the form 

v' (x) = T-'' + Jes X (x - x-'), xe D^, (2.4) 

where 63 is the unit vector perpendicular to the plane of motion. Therefore, if r = aei + 6e2, then 
63 X r = — 6ei + ae2. Equation (|2.4|1 shows that the velocity of a particle , j = 1,2, . . . , N is 
completely determined by two parameters: a constant translational velocity vector T-' and a scalar 
angular velocity uj^ . Both and to^ are unknown and must be determined in the course of solving 
the problem. 

Since each rigid disk is in equilibrium, the total force and torque exerted on by the fiuid 
must be zero, which provides the boundary conditions on the particle boundaries dD^: 

I Sn^ds = Oand [ x Sn^^'^ ds = 0, fov j = 1,2 ... N, , (2.5) 

where n-^ is the exterior unit normal to dD^ , and 

S = 2/ie(v) - PI. (2.6) 
In 1)2.6(1 and throughout the paper, e(v) denotes the strain rate tensor defined by 

e(v) = i(Vv + Vv^), (2.7) 

the superscript T stands for the transposed tensor, and / denotes the unit tensor. 

Solving equation (|2.1() with the boundary conditions l|2.2() (or <\2.'3\i ) and (|2.5|1 is equivalent to 
minimizing the functional 

over the function space U of admissible velocity fields u. This space is a space of vector functions 
in O satisfying either ()2.2|) or and such that 

n 

u = Bj, Uj G H'^ii^p), j = l...n, divu = 0, ((131) holds. (2.9) 

i=i 

Since the fiuid velocity v is the minimizer of the variational principle ()2.8|) - (|2.9() . the energy dissi- 
pation rate E in the fiuid (defined in equation (j3. 21 below) can be written as 

E = WnA^) = minl^f^^(u). (2.10) 

In next section we will see that calculation of the effective viscosities essentially amounts to 
calculation of E. 
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3 Effective shear and extensional viscosities. 



3.1 Effective dissipation rates 

We suppose that the suspension can be modelled on a macroscale by a single phase viscous fluid, 
called an effective fluid. The velocity field of the effective fluid is denoted by v° . The effective fluid 
is subject to the same external boundary conditions as the flow of the suspension. 

We assume that the effective stress tensor satisfies the constitutive equation of the form 

5° = F(e(vO)), (3.1) 

where F is a symmetric tensor function that does not depend explicitly on x. In this paper, we 
do not derive or postulate the precise from of the constitutive law for the effective fluid. Instead, 
we use the fundamental principle (going back to 8 ), that viscous energy dissipation rate of the 
suspension must be equal to the dissipation rate of the effective homogeneous fluid. The dissipation 
rates are defined by 

E= S ■ e(v)dx = 2fi e(v) • e(v)(ix (3.2) 

Jnp JQ.F 

in the suspension, and 

= f 5°-e(v°)dx, (3.3) 

in the effective fluid. In the equations ()3.2() . 1)3. 3() . S ■ e = SijBij is the inner product of tensors. 

For small particle volume fractions, ([H], H), this principle was further combined with the 
assumption that the effective fluid is Newtonian with a constant effective viscosity. However, for 
concentrated suspensions this assumption is not validated by rigorous mathematical derivation or 
experimental data and at present the problem of finding the effective constitutive law for such 
suspensions is still under investigation. Some experimental studies (e.g. JHl, ^ZI) suggest that 
the effective fluid may be non-Newtonian. Our calculations of the effective viscosity suggest that 
non-Newtonian behavior is possible for irregular (non-periodic or random) suspensions. 

We use the rheological definitions of shear and extensional viscosities as ratios of the corre- 
sponding components of the stress and strain rate tensors. To calculate the asymptotics of the 
two viscosities, we employ the network approximation introduced in (I I];). We analyze the network 
functional (discrete dissipation rate) introduced in ^ and show that the standard relation between 
two viscosities which holds for Newtonian fluids (see ^31 ch. 9 for the 3D-case and Appendix A for 
2D-case ) does not hold. 



3.2 Shear viscosity 

Suppose that a homogeneous effective fluid undergoes a steady shear flow with the shear rate 7. The 
velocity field v^^ satisfies the shear type boundary conditions 1)2. 2|) . The effective shear viscosity is 
defined by (see (|A.3|) ) 

= — = ^4^' ^3-4) 
7 7 i''i 

where is the corresponding component of the effective stress tensor and = dx. Note that 
our definition of in (|3.1|) implies that is constant when e(v'') is constant. Since E = F", the 
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equivalent definition is 

^i* = 2E^~'^\n\'^ = ifij-^\n\-^ [ e(v,0 • e(v,fe)dx (3.5) 
Thus the calculation of fi* amounts to evaluation of the total dissipation rate integral 

Esh = 2/i / e{vsh) ■ e{vsh)dx, (3.6) 

where Vgh solves (jTT |) -(P3 |) . 
3.3 Extensional viscosity 

A steady extensional flow of the effective fluid is characterized by a constant extension rate e. The 
velocity satisfies the extensional boundary conditions (|2.3|) . The extensional viscosity (see, e.g. 
[T^ . ch.9) may be defined by 

qO cO 

y = ^^n_^^ (3.7) 

where Sii, S22 are components of the effective stress tensor. Since = JqS^ ■ e(v)g^j(ix = {S^i — 
822)^1^1, the effective extensional viscosity can be defined in terms of the suspension dissipation 
rate E, as follows. 

E 



e^lf^l Jnp 

and calculation of A* again reduces to evaluation of the total dissipation rate (|3.6() with v^/i replaced 
by Vexi- In the remaining part of the paper we derive asymptotic formulas for the total dissipation 
rate under boundary conditions (|2.2|) and (|2.3j) . 



4 The Network 

Let us consider an arbitrary distribution of circular particles (disks) D^, whose centers are points 
X* in Q, for z = 1, 2, . . . , N. We suppose that is close to Amaxi so that neighboring particles 
can be close to touching one another. The network consists of vertices x* and edges. The edges 
connect only vertices that correspond to " neighboring " particles. Note that while for a periodic 
array the notion of a neighboring vertex (particle) is obvious, for non-periodic (e.g. random) arrays 
of particles it is not immediate. We introduce it via a Voronoi tessellation, which is a partition of a 
plane (or a planar domain) into the union of convex polygons Vi, called Voronoi cells, corresponding 
to the set of vertices x*. A Voronoi cell Vi consists of all points in the plane which are closer to x* 
than to any other vertex x.^ ,j 7^ i. 

The edges of Vi can lie either on d^l or in the interior of O. On each face of Vi, that lies inside 
0,, I X — X* 1 = 1 X — x-' I, for some i 7^ j. 

Definition 4.1 For each i = 1,2 . . . N , define the index set Mi by 

J\fi = {j G {1, 2 . . . , A^}, j 7^ i, such that Vi and Vj have a common edge} . 

We call a particle a neighbor of if j belongs to Mi that is the vertices x* and xP have a 
common edge in the Voronoi tessellation. 
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Note, that according to this definition, two particles are not neighbors if their Voronoi cells have a 
common vertex but do not share an edge. 

The minimal distance between neighboring particles and is given by 

(5ij =1 x^ -x^' I -2a. (4.1) 

We call 5ij inter-particle distances. If a disk is close to the external boundary 50, we define the 
particle-boundary minimal distance 5j by 

5' = d\si{^\dVt)-a. (4.2) 

To model the high concentration regime, we assume that 5ij and 5i satisfy 

5ij = 6dij, 6i = 6di, (4.3) 

where 5 is a small parameter in the problem, and the scaled inter-particle distances dij,di satisfy 

c < dij < 1, c<di<l, (4.4) 

with a fixed positive c independent of i,j. 

Definition 4.2 The centers x* of the particles are called the interior vertices of the network 
(graph) T, i = 1,2...N. The interior edges bij of the network connect neighboring vertices x* 
and x-^. If a Voronoi cell Vi has edges belonging to dVf^yjdVt' , the corresponding x* is called a 
near-boundary vertex. Each near-boundary vertex x* is connected with dO.^ by an exterior 
edge, which is a segment hi perpendicular to dO,^ . This segment intersects dO.^ at an exterior 
vertex denoted by x*^*^ . 

Finally, the network (graph) T is the collection of all interior vertices, exterior vertices, and 
all the edges connecting these vertices. 

The set of indices of the near-boundary vertices will be denoted by /, that is i G / if the vertex 
x' is connected to dVi. Also, we write i G if x* is connected to dQ~^ {dil.~). 

Note that T is essentially the Delaunay graph dual to the Voronoi tessellation, and the above 
notions admit straightforward generalization to three dimensions. 

5 Network approximation of the effective viscosity 
5.1 Network equations 

To define the network approximation, we first assign a translational velocity T* and an angular 
velocity of a particle to the corresponding interior vertex x*. At the exterior vertices x* we 
prescribe the velocity vector g which represents the boundary conditions 1)2.2(1 or (|2.3|) . 

For each pair of neighboring particles and we introduce a gap Hij which represents a 
fluid region where lubrication effects are very strong as shown on Fig.l. The width Rij of such 
gap is independent of 5. For technical reasons it is convenient to work with non-intersecting gaps. 
Since the maximal number of gaps adjacent to each particle is equal to the maximal coordination 
number (number of neighbors), K, we can choose Rij = KM a, where M is sufficiently small and 
fixed. The precise value of M is not essential for our purposes. 
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Figure 1: A gap Iljj between the neighboring particles and DK 



T* 




Figure 2: Assignment of T*,u;* and orientation of the gap between two neighboring particles. 



The orientation of each interior gap Wij relative to a disk is specified by a unit vector 



|x* — 

We also let p*-' be the unit vector obtained by rotating q*-' clockwise by 7r/2 (see Fig.2). 

Boundary edges and corresponding gaps are oriented perpendicular to ei. This reflects the 
physical fact that the zone of the largest energy dissipation is located near the shortest line con- 
necting X* with the boundary. Therefore, q* = 82, (respectively —62), when I?' is adjacent to 
(respectively d^~). 

Next, to each edge of the network we associate a dissipation rate W^^ [W"^), calculated in the 
corresponding gap 11^ (Ilj). The calculation of the dissipation rates employs lubrication approx- 
imation in the gap. The velocity in the gap is decomposed into three velocities, representing the 
"elementary" motions called spring motion, shear, and rotation (see Fig. 3). The total velocity field 
in a gap is the sum of these elementary velocities and a "residual" velocity field, whose contribution 
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(5.1) 



spring motion shear motion rotation 




Figure 3: Three elementary motions. Arrows represent the boundary conditions. 



to the gap dissipation rate is 0(1) as 5 0. Lubrication approximations for each of the elementary 
velocities and estimates for the residual are obtained in j3j where more details can be found. 

Using approximations of the elementary velocities to calculate (up to the terms of order 0(1) 
as (5 — > 0) the dissipation rates in each gap we obtain 

^ 5-3/2^7]^ [(T^-Ti)-q^J]' + 

in the interior gaps Iljj, and 

= 5--"^Cl^ [(T* - g) • qf + r^/2c;,,a2 {2u;f + r^/^C:, [(T^ - g) • + au/f (5.3) 

in the boundary gaps Ilj. The expressions for factors Clp,Clj^ and C*^^ and Clp, Clf^, C^^^ are 
calculated explicitly in j4j: 

r<ii _ 9„,,/a\V2 _ 9„,,/'aNl/2 

Wot - TE^^^\.d~! ■ ^rot - TE^Hdl) ■ 

In the formulas 1)5. 4() , dij , di are the scaled inter-particle distances defined in (|4.4() . 
The sum of the local dissipation rates W^^ , is a quadratic form 



Q = En., W^^^' + En,: 

N 

i=l j e AT; 

+ ^''^^Ci^ [(T^-g)-p^ + au;f }, 
where I denotes the set of indices of the near-boundary vertices defined in Section 4. 



(5.5) 
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The form Q is the discrete network approximation of the functional W^j, in the variational 
principle (|2.8j) . The main idea of the network approximation is that the most of the energy is 
dissipated in the gaps Iljj , Ilj , so that 

E = Enet + 0{l), as (5^0. (5.6) 

The discrete dissipation rate Enet in (|5.6|) is defined by 

Enet = min Q = (5(T^j^, T^j„, ...T^j^, ...w^j^), (5-7) 

where the vectors T^^,^, i = 1, 2, . . . , and scalars u)"^ ,i = 1,2, ... ,N minimize Q in (|5.5)) . The 
minimum in ()5.7|) is taken over all possible collections of T* , . 

It is well known that solving the minimization problem for Q is equivalent to solving the linear 
system (Euler-Lagrange equations), which is obtained by equating the gradient of Q to zero. Setting 
to zero partial derivatives with respect to T^*, Z = 1, 2 we obtain 

J2 {'^"^^^Csi [(T* - T^') • P*^' + acu* + aJ] p*^'} + = F*, (5.8) 

for each i = 1,2, N, where 

= I '5-3/2cyT^•q^)q'+rV2cj,[T^.p^ + au;']p^ if i G / 
\ otherwise, 

(5.10) 

\ otherwise, 

(5.12) 

Next, equating the partial derivatives to zero we obtain 

E [(T' - T^) • p*^- + ao.' + au:^] } + 

E {s-'/'CZt i-^' - }+l3' = M\ (5.13) 
for all i = 1, . . . , A^, where 

1 otherwise. 



M = \ f/'f-(g-P^) if^^^' (5.15) 
' otherwise. 
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Equations (|5.8() and (|5.13l) are, respectively, the equations of force and torque balance of the 
particles, and the minimization in (|5.5j) ensures that the rigid body translational and angular 
velocities are chosen in such a way that the suspension is in mechanical equilibrium. Note also 
that (|5.8j) is a system of 2N equations, and (|5.1,3j) is a system of equations. Together they form 
3A^ equations for 3A^ unknowns (T*,a;'). The coefficients and right hand side of (|5.8() are of order 
(5~^/^ and 5"^/^ while all the coefficients in (j5.1.S|) are of order When all T* are zero (no 

translations), the remaining terms are of order (5~^/^, but in the case = 0,i = 1,2,... (no 
rotations), the remaining equations contain terms of order (5~^/^. This reflects the well known fact 
that the contributions from local translational spring motions are stronger than the contributions 
from rotational and other translational motions. Keeping only spring translations we obtain the 
truncated (leading order) discrete dissipation functional 



1 ^ 



(5.16) 



=1 j&Ni i&I 

Then Q can be decomposed as follows. 

Q = S-^/^Q + (5.17) 

where the coefficients of the forms Q and Q' do not depend on 6. 

We next show that the total discrete dissipation rate Enet can be estimated by the truncated 
dissipation rate obtained by minimizing Q. Let T', i = 1,2, . . . , N denote the translational velocities 
which minimize Q ( they are clearly independent of 6). Since Q' in (|5.17|) is non-negative. 

Since Q'(t\u:^ = 0) < C with C independent of (5 as (5 ^ 0, (OKI) implies 

Enet = &''^^'^E + 0{5'^'^), as 5^0, (5.19) 

where 

E = Q(X)= min g(T\...T^). (5.20) 

This equation enables one to calculate the leading term in the asymptotics of the effective viscosity 
by solving a simplified minimization problem involving only the translational particle velocities. 
However, this algorithm is useful only when 

^ > 0, (5.21) 

because in this case the leading term in the asymptotics of the dissipation rate is of order 5~^/'^ . If 
minQ = 0, the leading term degenerates, and the rate of blow up in (|5.2()|) is at most (5~^/^. 

Minimization of the truncated quadratic form Q corresponds to solving the truncated linear 
system of the Euler-Lagrange equations 

[C%{T^' - T^') • q^^q*^] + B(? ) = R^ i = 1, 2, . . . , iV, (5.22) 
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where 



Cjp(T» • q^)q*, when i G / 
otherwise. 



B(T') 

The right hand side vectors R' represent external boundary conditions: 



(5.23) 



1 otherwise. 



(5.24) 



To better see the structure of the functional Q and the system (|5.22j) it is convenient to rewrite 
them in a compact form as 

Q(z) = -Az-z-f -z + r, (5.25) 



and 



Az = f, 



(5.26) 



where z is the particle velocity vector, and f is the vector of discretized boundary conditions. The 
vectors z, f G R^^ are defined by 



^2 

^2 



V Ti^ J 



( ^\ \ 

2 



and 



^1 

i?2 



(5.27) 



The fixed scalar r equals f • f where components of f are defined by the components of f , as follows. 



2k-l 



f2k- 



f2k 



f2k k 



1,2,, 



.N. 



The matrix A in ()5.25|) is symmetric and non-negative definite. Its entries are determined by the 
scaled inter-particle distances dij , particle radius a and vectors q'-' , q' defined by 1)5. 1|) . The vector 
f and the scalar r in ()5.25|) are determined by the boundary conditions on dQ^,dQ~ (extensional 
or shear). The linear term f • z depends on the translational velocities of the particles located near 

dn+ or dn-. 

To the leading order in 6, the effective viscosity is determined by the discrete dissipation rate 
^-3/2^ where E is the minimum of Q. Thus the qualitative behaviour of the effective viscosity is 
determined by solving ()5.26() . We now sketch the issues which arise in computing the leading term 
of the effective viscosity. First, in the case of shear viscosity, f = 0, which means that the shear 
boundary conditions do not contribute to the strong blow up. This results in the weak blow up of fi* 
(see section 6). In case of extensional conditions, f 7^ 0, which means that the right hand side of the 
full network equations 1)5. 8|) contains terms of order (5~^/^. In this case, calculation of the leading 
term is impeded by the fact that the matrix A is not invertible. Indeed, from (|5.16)) it is clear 
that the value of Q does not change when all vectors T*,i = 1,2,... are replaced by T* + tei 
(horizontal translation by tei, t arbitrary real). This is not surprising, because the functional Q 
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is invariant under horizontal translation of T*. The invariance is due to the vertical orientation 
of the boundary gaps explained in section 5.1. Note that this is not equivalent to translation of a 
coordinate system, since g is not changed. Translational invariance of Q implies that any z* that 
solves the non- homogeneous system (|5.26|) produces other solutions of the form z* + iwo, where t 
is arbitrary real, and 

wo = (1,0, 1,0,..., 1,0)^. (5.28) 
This means that vectors of the form iwg solve the homogeneous system 

ylz = 0. (5.29) 

Because of the above mentioned invariance of the functional Q, the leading term in the asymptotics 
of effective viscosity will be uniquely defined, unless (|5.29|) has some other nontrivial solutions. 
Therefore, it makes sense to look for conditions on the network which would guarantee that every 
solution of the homogeneous system is of the form twQ. Then the non- homogeneous system (|5.26|) 
would be uniquely solvable up to horizontal translation. In Section 7.3 we show that for a typical 
random distribution of particles this is indeed the case. 

The last issue concerns the validity of the estimate (|5.21j) . The functional Q is non- negative, 
but it may be zero. When (|5.2H) holds, local lubrication analysis provides the correct order of the 
leading term in the asymptotics of the extensional effective viscosity (5~^/^ in dimension two and 
5^^ in dimension three). If ()5.21|) does not hold, that is, 

minQ = 0, (5.30) 

then the leading term in dimension two is of order 5~^/'^ , (ln(l/5) in dimension three). 

Whether or not the estimate 1)5. 21() holds, depends on the geometry of the particle array as well 
as the boundary conditions on and dQ.~ . In section 7.4 we show that in the case of extensional 
boundary conditions, (|5.21() holds for generic arrays, and thus the leading term in the asymptotic 
of Enet ( and extensional effective viscosity A*, see (|3.8|) ) is of the order However, there exist 

special arrays for which the extensional effective viscosity is of order 5~^/'^ . An example of such an 
array is presented in section 7.5. 

6 Effective shear viscosity 

In this Section we show that in dimension two, the asymptotic order of the shear effective viscosity 
/i* is while the local lubrication analysis predicts the rate 5~^/'^. In three dimensions, 5~^/'^ 

and 5^^/'^ should be replaced by, respectively, In 5 and 5^^ (see ^). The local analysis in three 
dimensions predicts that the asymptotics of the shear effective viscosity (see ()3.5() 1 should be 
of order , but numerical simulations in and experimental results in |15j and |12j show that 
random suspensions in shear flow have effective viscosity of order In (5. Our estimate ^* = 0(5"^/^) 
is therefore in agreement with the three-dimensional results in up to the difference in the 
dimension of the space. 

The decrease in the asymptotic order of /x* is a global phenomenon, which shows that the 
local analysis could be misleading, and analysis of the entire particle array leads to qualitatively 
different results. This fact can be explained as follows. The "strong" blow up rate ((5~^ in three 
dimensions and 5^^^"^ in two dimensions) is obtained using classical lubrication techniques applied 
to two particles D^,D^ whose translational and angular velocities T', T-', w', a;-' are independent 
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Figure 4: A two-disk network with shear boundary conditions. 



of 6. However, the network analysis of the ensemble of particles, interacting with each other and 
with the external boundary, shows that the above velocities may depend on 6. In the case of 
shear boundary conditions, the right hand side of the network equations (|5.8() . H5.13() is of order 
whereas the matrix of the network equations is of order 6~^^'^. The solution of the network 
equations is thus small (of order 6) which makes local and global dissipation rates small. 

When the boundary conditions are given by (1221), the vectors R* in KWi are zero (since 
q* • ei = for all i = 1,2,... ,A^). Consequently, the right hand side f in ()5.26|) is zero. The 
functional Q reduces to Az ■ z which is clearly zero for every solution of the homogeneous system 
Az = 0. Since z = is an admissible trial vector for Q, E = min Q = 0, and from (|5.19j) we obtain 

1^' = CEnet = 0{6~'/^), (6.1) 

where C = 2-f~'^\n\-^ (see ((HSl). The same conclusion can be obtained by directly estimating 
the minimum of the form Q. Consider first the shear boundary conditions (|2.2I) with 7 = 1 and a 
simple example of a two-disk network on Fig. 4. The functional Q for this example has the form 

Q =5-'" [cll [(Ti - T2) . qi2]2 + [(Ti - ei) -62]^ + 0% [(T^ + e^) • e^]^} 

+^-'/' [(T^ - ei) • ei + a^^f^ [(T^ + e^) • ei + aa;^]'} 

The dissipation rate E is the minimum of Q. Hence, for any collection T^, T^, wi, 0^2, we have 
E < (5(T\t2,u;i,u;2)- In particular, choosing = = 0, wi = u;2 = in we obtain 

E<6-'/\Cl, + C!,) 

Since Cj^, C^^ are independent of 6, the blow up rate of E is at most 5^^/^. 
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Next, consider the general case. Since for shear boundary condition g-Lq* for all i G /, we 
obtain 

N 

1=1 j G AAi 
j < i 

N 

i=i . GAA, ^ ^ (6-3) 

Choosing the trial vectors T* as T* = 0, i = 1, 2, . . . , and = 0, i = 1, 2, . . . , we obtain 
Er^et = min g < Q(T* = 0, o;* = 0) = r ^ ^^.(g • p^)^ = 72^-1/2 ^ c*:,. (6.4) 

where we used the shear boundary conditions ()2.2|) . The sum in the right hand side of 1)6.4(1 is 
independent of 5 and the shear rate 7. This estimate shows that under the boundary conditions 
(|2.2|) . the effective viscosity blows up as 5'^/"^ at most. More precisely, we have the following 
Conclusion. Effective shear viscosity of a concentrated suspension has the blow up rate 5^^/^ in 
two dimensions, that is 

< Cb-^l'^ as 5 ^ 0, 

with C independent of6,'y. 



7 Extensional effective viscosity 

7.1 Simplification of boundary conditions and outline of the method 

In this Section we show that for the extensional boundary conditions, the leading term in the 
asymptotics of the dissipation rate E from (|5.19() may or may not be zero depending on the geometry 
of a particle array, that is, the leading term in the asymptotics of the effective extensional viscosity 
is either of order 5~^^'^ (strong blow up) or 5~^^'^ (weak blow up). We provide two geometric 
conditions on the network graph which ensure strong blow up. 

In a planar steady extensional flow of the effective fluid, the rate of strain tensor is 

e(v») = (j (") 

where e denotes a constant extension rate. The corresponding velocity field is of the form 

v0 = (exi,-ex2)^, (7.2) 

which gives the boundary conditions 

^0 _ I (ea;i, -e)"^, when X2 
[(exi,e)"^, when 2:2 = 



1 (on9n+), 

-1 {on on-). 



(7.3) 
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We decompose v*^ as 

vc ' ^ ge 

where vj]^ is a vertical contraction velocity satisfying 



v° = vS, + vO, (7.4) 



j-ee2 on 80.+ , 

Kc = g^>c = < (7.5) 
I ee2 on asZ , 

and v^g is the horizontal extension velocity field with the boundary conditions given by 



Since 



Jea;iei on 80.+ , 

'ge = Sge < (7.6) 

exiei on aiZ . 



g = gi-c + gge, (7.7) 



and gge-Lq*,^ G I, the value of Q in ()5.16|) does not change when g in ()5.16|) is replaced by gyc- 
Hence, 

Q{T\io\g) = 5-3/2Q(T^ g,e) + 6-'/^Q'{T\u;\ g,, + g^^), (7.8) 

To determine the rate of blow up of the dissipation rate, we need to analyze the minimizers of 
Q(T',g^c)- The estimate (|5.19() implies that the second term in the right hand side of H7.8() is of 
order (5~^/^ (at most). Since Q is independent of 6, its minimizing vectors T* are also 5- independent. 
Consequently, the blow up rate of the dissipation depends on whether the minimum of (5(T*,gt,c) 
is positive. If it is, the extensional effective viscosity A* is of order otherwise A* grows no 

faster than Which type of behavior occurs, depends on the validity of the estimate (|5.21|) . 

As mentioned in section 5, (|5.2H) may fail for certain particle arrays. In this section we provide 
geometric conditions which insure positivity of minQ, and give examples of networks for which this 
minimum is zero. The principal conclusion here is that extensional viscosity of suspensions with 
a comparable volume fraction of particles may vary by an order of magnitude in the inter-particle 
distance, depending on the geometry of a particle array. 

Our method of analysis is based on the following simple observation. The form Q(T*,g^c) is a 
sum of non- negative terms, namely 

1 ^ 

Q(T^ g,,) = 2 E E - • q'')' + E ^spiiT^' - g-) • q')'- (^.9) 

i=l j&Afi i€l 

This shows that min Q(T*, gt,c) = if and only if the minimizing vectors ,i = 1,2, . . . , N satisfy 
the system of equations 



(T^ - T^) • ci'^ = 0, i = l, 2..N,j G 
(T* - g^e) • = 0, 



(7.10) 



Hence, if ()7.1U() does not have solutions, (|5.21l) must hold. Also, it should be noted that if the 
estimate (|5.21|) holds for some subgraph of the network, then it also holds for the whole network. 
This observation can be used to reduce the network to a simpler graph, for which it is easier to 



prove (|5.2H) . 



16 




I • 1 

Figure 5: A network of four vertices with E = 0. 



7.2 Simple examples 

Suppose that the boundary conditions are given by H7.5j) with e = — 1. In this section we present 
three simple examples of networks with small number of vertices. One of these networks has E = 0, 
and for the other two examples E > 0. 

Example 1. This an example of the network for which E = 0. To demonstrate this, we show that 
there is a nontrivial particle velocity vector z such that Q(z) = 0. Consider the network on Fig. 5. 
The vectors q*-' are defined as follows. 

qi = 62, q^ = -62, 

q^^ = q^^ = 75(ei -62), (7.11) 

Next, define T* as follows. 

T^ = e2, T2 = -ei, = ei, = -62. (7.12) 
The functional Q corresponding to this network is 

Q = Ci,[(Ti-e2).qi]' + Ci2[(Ti-T2).qi2]2 + 

[(Ti - T3) . qi3] 2 + _ T^) . q24] 2 + (7.13) 

CU [(T3 - T^) . q34] 2 + c% [(T^ + 62) . q4] ^ 

When T* are defined by (|7.12|) . all the scalar product in brackets in (|7.13|) are zero, and therefore 
minQ = 0. 

Example 2. For the network of three vertices in Fig. 6, minQ > 0. To show this, consider the 
system corresponding to the general system 1)7. 10() . 

(Ti+62)-e2 = 0, (Ti - T2) . qi2 = 0, (T^ - T^) • qi3 = 0, 

(T2 _ T3) . q23 = 0, (T2 - e2) • 62 = 0, (T3- 62) -62 = 0. ^'''^^ 

We prove that the system (|7.14j) has no solutions. Indeed, the last two equations imply = 
t26i + 62, T3 = ^361 + 62, for some scalars t2, t^. Next, the first equation in the second row of 1)7. 14(1 
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Figure 6: A network of three vertices with E > 0. 



Figure 7: A network of four vertices with E > 0. 



yields {t2 - h)ei ■ q^^ = 0, and thus t2 = ts = t. Substituting 
and third equations in the first row of H7.14() we obtain 



(Ti - tei 



e2j • q 

62) 



12 

,13 



0, 
0. 



tei + e2 into the second 



(7.15) 



Since q^^,q^^ are non-colhnear, (|7.15j) yields = tei +62, which contradicts the first equation in 
the first row of H7.14() . 

Example 3. Next, consider a rectangular network of four vertices in Fig. 7. The system (|7.10|) 
for this example becomes 



(Ti - T^) . ei 
(Ti +e2) -62 : 



_ 0^ - T^) • e2 = 0, 

0, (T2 + ea) • 62 = 0, 



(T2, 
(T3 



T3) . e2 
62) • e2 -- 



0, 



(T3 



-n-ei 

62) • 62 = 



= 0, 
0. 



(7.16) 



The last two equations in the second row of (ITTKl) yield = tsei + 62, T'' = Uei + 62, with some 
scalars ta, ^4. Next, the second and third equations in the first row of (|7.16j) produce = tiei +62, 
= t2ei + e2, which contradict, respectively, the first and second equations in the second row. 
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Figure 8: An example of the network. The shaded quadrilateral represents a defect cell. 



The three examples above seem to indicate that two basic building blocks for networks with 
E > (strong blow up) are triangles (Fig. 6) and rectangles aligned with the edges of (Fig. 7). 
Misaligned rectangular structures such as shown in Fig. 8 would produce E = (weak blow up). 

7.3 Quasi-triangulated graphs 

7.3.1 Definition of quasi-triangulated graphs. Solvability of the system ()5.26|) 

The network F partitions 0, into a a disjoint union of convex polygons, which are called Delaunay 
cells. When points x-' are distributed randomly in 0, the interior Delaunay cells are typically 
triangles. This simple but important fact can be explained as follows. The edges of Voronoi 
tessellation are perpendicular bisectors of the edges of Delaunay cells. If an interior Delaunay cell 
is, for instance, a quadrilateral, then any two vertices lying on a diagonal cannot be neighbors, and 
therefore the point of intersection of four edges of the Voronoi tessellation must be equidistant from 
the four vertices. This means that a convex quadrilateral may be a Delaunay cell only if all four 
vertices lie on a circle. When the vertices of the network are randomly placed, the likelihood of four 
(or more) points lying on the same circle is small. It is natural to call such cells the defect cells. 
Therefore, most of the interior Delaunay cells of a random network are triangles. Other polygonal 
cells (quadrilateral, pentagonal etc.) are typically small in number, isolated and are likely to be 
unstable in the actual flow. 

Since the exterior edges of the network are vertical, the cells adjacent to the boundary are 
typically quadrilateral. An example of a generic network is shown on Fig. 8. 

Next, we define a broad class of graphs containing generic Delaunay- type networks. The graphs 
in this class are called quasi-triangulated, because of the presence of triangular cells. The number 
of triangular cells may be relatively small, as indicated by the examples below. We show that 
for quasi-triangulated graphs the leading term in the asymptotics can be effectively computed by 
minimizing the functional Q, that is the minimum is unique and can be found by solving the linear 
system (|5.2()|) . We also show that the crucial estimate (|5.21|) holds even for more general graphs 
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t t t t t 

Figure 9: A rectangular graph on the left is not quasi-triangulated. Clearly, in this case Tp = T~ . 
The triangular graph on the right is quasi-triangulated. 



which contain a "spanning" quasi-triangulated subgraph. 

Given an arbitrary network graph T, we define its maximal quasi-triangulated subgraph Tm by 
the following iterative procedure. 

Step 1. Consider interior vertices which are connected to dVt~ and call these vertices generation one 
vertices. All interior edges connecting these vertices are generation one edges. Add all generation 
one edges and vertices to the subgraph. 

Step 2. Consider all remaining vertices which are connected to the vertices of the subgraph by at 
least two non-collinear edges. These vertices and edges belong to generation two. Add them to the 
subgraph. Note that the non-collinearity condition leads to formation of "supportive triangles". 
Step 3. Repeat step 2 until no more vertices can be added. 

If the maximal quasi-triangulated subgraph Tm contains all interior vertices of F, we call the 
graph F quasi-triangulated. It turns out that for quasi-triangulated graphs the system (|5.26|) is 
"almost uniquely" solvable, that is two solutions differ by a vector of the form two, where wq is 
defined by (|5.28|) . and the value of the form Q for these solutions is the same. 

Proposition 7.1 Suppose that the network graph F is quasi-triangulated. Then there is a unique 
solution of the system h5.2b}) . up to a horizontal translation. 

This proposition is proved in Appendix B. 
7.3.2 Examples of quasi-triangulated graphs 

First we observe that a restriction to of a periodic rectangular lattice is not quasi-triangulated. 
By contrast, a periodic triangular lattice restricted to Q. is quasi-triangulated, (see Fig. 9). 

In fact, if a network F is not periodic, but all of its interior Delaunay cells are triangles, then F 
is quasi-triangulated. The converse is false (see an example in Fig. 11). This means that the quasi- 
triangulation property is more general than triangulation property. In Q we considered a class of 
graphs containing a triangulated path (see Fig. 10). In Fig 11, we present two examples showing 
that a quasi-triangulated graph need not contain a triangulated path. On the other hand, a graph 
containing a triangulated path is not necessarily quasi-triangulated. However, when a triangulated 




20 



Figure 10: A triangulated path. 




path exists, it must be contained in the maximal quasi-triangulated subgraph. Therefore, in this 
case the maximal quasi-triangulated subgraph is spanning ("extends from top to bottom"). Below 
in section 7.4.1 we show that this condition is sufficient for positivity of the minimum of Q. 

7.4 Strong blow up of A*. Percolating rigidity networks 
7.4.1 Quasi-triangulated subgraphs 

For the rest of this section, we consider the steady flow of the suspension corresponding to the 
boundary conditions 1)7. 5|) with the extension rate e = 1. 

We shall say that a network is a percolating rigidity graph when the crucial estimate (|5.21|) 
holds. In this subsection we show that quasi-triangulated graphs are percolating rigidity graphs. 
Moreover, a graph has percolating rigidity even when it is not quasi-triangulated but contains a 
spanning quasi-triangulated subgraph. 
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Proposition 7.2 Suppose that the boundary conditions are given by i7.5\ ) and the network graph 
r contains a spanning quasi-triangulated subgraph. Then T is a percolating rigidity graph. Con- 
sequently, the extensional effective viscosity \* of suspensions corresponding to such networks is 

The proof of this proposition is given in Appendix C. 

7.4.2 Networks containing a vertical path. Periodic square networks 

We now present another class of percolating rigidity graphs. Namely, these are graphs that contain 
a path connecting dVt^ and dVt^ ^ such that all edges in this path are oriented along e2-direction 
(vertical). The simplest representative of this class of graphs is a periodic square lattice. A periodic 
square graph does not contain a spanning quasi-triangulated subgraph (in this case Tm is just a 
path which consists of the vertices adjacent to dVt~ , see Fig. 9). However, the estimate (|5.21j) holds 
for a periodic square graph such as the one shown on Fig. 9. This follows from the more general 
criterion 

Proposition 7.3 Suppose that a network graph F contains a path F^ such that 

i) it connects and dVt~ , 

ii) all edges ofV^ are vertical. 

Then T is a percolating rigidity graph. 

As noted above, if a path F^ has percolating rigidity, then the "larger" graph F is also a 
percolating rigidity graph. Therefore, we only need to show that a path of vertical edges has 
percolating rigidity, that is, the quadratic form Q corresponding to such path is positive-definite. 

For simplicity, consider a path containing three vertices. The argument can be directly gener- 
alized to an arbitrary number of vertices. As before, ()5.21|) will hold if the corresponding system 
(|7.1()|) has no solutions. The latter now has the form 

(Ti-e2)-e2 
(Ti - t2) . e2 
(T2 _ t3) . e2 
(T3 + es) • 62 

Introduce new unknown vectors T* = t*+e2 i = 1, 2, 3. The last equation of (|7.17l) yields = tsei, 
where t^ is a scalar. Then from the third equation we obtain = ^2^1, and then the second 
equation yields = tiei, which contradicts the first equation. The contradiction shows that the 
system 1)7. IZp has no solutions. The argument can be easily modified to show that if at least one 
of the edges of a path is non-vertical, then this path does not have percolating rigidity. 

Proposition 17.31 shows that for rectangular lattices oriented parallel to the edges of dQ, asymp- 
totics of extensional effective viscosity is of order 5~^^'^ which corresponds to strong blow up. 
Asymptotic formulas in also predict strong blow up of the extensional viscosity for cubic lat- 
tices. Thus our results are consistent with the results in We refer to section 8 for a more 
detailed comparison. 



0, 
0, 



(7.17) 
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Figure 12: A rotated rectangular lattice of 12 vertices. 



7.5 Weak blow up 

In this section we present an example of a particle array for which the leading term in the asymp- 
totics of A* is zero. Roughly speaking, the array in question is a rectangular lattice rotated so that 
none of its interior edges are vertical. Let k denote a unit vector non-collinear to either ei or e2. 
The interior edges of the rectangular network in Fig. 12 are either parallel or perpendicular to k, 
while the prescribed boundary velocities are parallel to e2. This misalignment will lead to weak 
blow up. To show that the system 1)7. 10() is solvable, we first consider a single path connecting dQ~^ 
and dQ^ . This path may be any of the three such paths in Fig. 12. The argument we use admits 
a straightforward generalization to a network with an arbitrary number of vertices. 
The system (|7.10|) written for the path has the form 



(7.18) 



For technical reasons it is convenient to introduce new unknowns ui = — e2, ui2 = — T^, 
U23 = - T^, U34 = - T^. The relations between u^^ and T* are 

= ui + 2e2, 
T2 = ui + 2e2-ui2, 
= ui + 2e2 - U12 - U23, 
T'' = ui + 2e2 - U12 - U23 - U34. (7.19) 

From the first four equations of ()7.18() we see that ui = tiei, and Uj^j+i = tj^j+ik-*-, i = 1,2,3, 
where k"*" denotes a unit vector orthogonal to k, and are scalars. The problems of solving 

(|7.18() is now reduced to finding ti, ij^j+i such that the fifth equation of (|7.18|) is satisfied. The fifth 
equation shows that 

= t4ei - 62, (7.20) 



(Ti 


- 62) 


62 = 


0, 


(Ti 


_T2) 


•k = 


0, 


(T2 


_T3) 


•k = 


0, 


(T3 


_T4) 


•k = 


0, 


(T^ 


+ 62) 


62 = 


0, 
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where is a scalar. Equating (|7.2UI) and (|7.19() we obtain the equation for ti, ti2, ^23; ^34 and 

Uei = tiei + 62- (ti2 + t23 + t34)k^. (7.21) 
This yields two scalar equations 

= l-(ti2 + t23 + t34)k^-e2, (7.22) 

and 

ti-t4 = (ii2 + t23 + i34)k^-ei. (7.23) 

The system of two equations (|7.22() . H7.23|) for five unknowns has infinitely many nontrivial solutions 
as long as k"*" • e2 7^ , that is, the interior edges of the path are non-vertical. 

At the next step of construction, we consider the full lattice on Fig. 6 containing 12 vertices. 
The interior edges of the graph are oriented either by the unit vector k as above (longitudinal edges), 
or by k"*- (latitudinal edges). We view the graph as the union of three paths of longitudinal edges 
extending from dfl~ to dft^, with latitudinal edges connecting these paths. To obtain the desired 
example, we choose the vectors T* for one of the paths, as explained above. Then we prescribe the 
same vectors to the corresponding vertices of two remaining paths. Now, if two neighbors x* , x-^ 
belong to different paths, then T* = T-' , so the corresponding equation (T* — T-^ ) • k-*- = of the 
system (|7.10|) is satisfied. When the neighbors x^jX-' belong to the same path, the corresponding 
equation (T* — T-') • k = is satisfied by the choice of T*,T-'. Therefore, the whole system 1)7. 10() 
in this case has infinitely many nontrivial solutions. Each of these solutions makes Q, and thus the 
leading term in the asymptotics of A*, zero. 



8 Comparison with some results for periodic cubic arrays in di- 
mension three 

The main objective of network approximation is to define effective properties for non-periodic 
deterministic or random arrays. While techniques of periodic homogenization are well developed 
(0) ^2j; ^Oj and references therein), non-periodic geometries are much less understood. In the 
end of this section, we compare our results applied in particular case of a periodic square array (in 
dimension two) with the results of JJ| obtained for cubic arrays in dimension three. 

The effective viscosity of an infinite periodic suspension obtained in 11 is the fourth order tensor 
fj,* (in this section we use the notation from from TT ) . In an effective flow with the constant strain 
rate 7, the effective stress is 

''■'ij — ^H-ijkl iki - ^ ^ij, 

where P is an effective pressure. The authors of 11 obtained the following formula for the com- 
ponents of n* : 

12 1 
t^tjki = 2^(1 + P)i^ikSji + 5u5jk - -6ij6ki) + fj.{a - I3){5ijki - -Sijdki), (8.2) 

where Sijki = 1 if all indices are equal, otherwise Sijki = 0. ;U is the fluid viscosity, and a, (3 are 
functions of the small parameter e, related to the inter-particle distance 6 as follows: 

e = ^.A, ,,,^0. (8.3) 



50. = 2^^,,,7fc,-P5,,, (8.1) 
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When the effective fiow is incompressible, = 0, so the formula (|8.2() simplifies to 

1 

f^ijki = 2/^(1 + (^)i^ik^jl + ^il^jk) + ^(a - P)Sijki. (8.4) 
For simple cubic lattices, up to the terms of order 0(1) in e, (XL) 



a = ^7re i + ^vrlne \ (8.5) 

and 

/3=-7rlne~\ (8.6) 

Suppose that the imposed effective flow is a steady shear with the velocity v = {kx3,0,0) 
(three-dimensional analogue of HA.1|) ). where k > is a constant shear rate. The components of 
the corresponding strain rate tensor 7 are 713 = 731 = k > and ^ij = for other values of (^, j). 
Then, using 1)8. 4() . we obtain from ()8.1() that the only nonzero components of the effective deviatoric 
stress 2/2*7 

(2/i*7)i3 = (2^^*7)31 = 2fi{l + P)k. (8.7) 

Since a is not present in (|8.7)) . the nonzero components of the deviatoric effective stress are of 
order lne~^ and thus the shear effective viscosity calculated by the three-dimensional analogue of 
definition (|.S.5|1 is of order lne~^ (weak blow up in dimension three). The weak blow up was not 
identified in ^TTj, but it can be easily deduced from the formulas derived there. 

In the case of an extensional flow, velocity vector v = {kxi, kx2, — 2^x3), (compare with HA.5|l ). 
where k > is a constant extension rate. The strain rate tensor is 




7 = K (8.8) 
\ -2k J 

Then, using (|8.1|) and (|8.4j) we obtain the deviatoric stress 

/ L 

2/x*7= L I (8.9) 





where 

L = 2Kfi{l + [5) + Kn{a- (3). (8.10) 

Components of the deviatoric stress contain a and are therefore of order e~^. Consequently, the 
extensional effective viscosity is of order (strong blow up in dimension three). 

Although Noonan and Keller did not address the issue of weak versus strong blow up for the 
effective viscosity in the formulas ()8.5|) - 1)8. 1U() are consistent with the results for square lattices 
in Section 6 of this paper in the following sense. If a periodicity cell corresponding to a simple 
cubic lattice is subjected to a uniform shear (extensional) flow, then the straightforward calculation 
presented above shows that the shear (extensional) effective viscosity exhibits weak (strong) blow 
up. However, for other lattice types such as BCC and FCC, formulas from jJJ imply strong blow 
up of both viscosities, while our approach leads to the weak blow up of the shear viscosity for all 
two-dimensional lattices. Perhaps, this can be attributed to the fact that in effective viscosity 
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is defined in a different way, as explained in the introduction. In particular, our definition takes 
into account boundary effects which are known to be essential in rheological measurements (||17j. 
jl5j). Considerations based on infinite periodic lattices cannot capture these boundary effects. 
Futhermore, if our approach is applied to a rectangular periodic lattice in a rectangular domain 
chosen so that lattice periodicity is preserved, then our results are consistent with 

9 Conclusions 

We have studied properties of the network approximation of the shear effective viscosity ^u* and 
extensional effective viscosity A* of two-dimensional flows of concentrated suspensions with complex 
geometry. The flow domain was chosen to be a square, upper and lower sides of which represented 
the physical flow boundary, while the other two sides were assumed to be free surfaces. 

A small inter-particle distance parameter 6 was used to describe the high concentration regime 
for particle arrays which are not necessarily periodic (i. g. random). In the recent paper ^1] by 
Sierou and Brady, the high frequency dynamic viscosity of concentrated suspensions (see ^21) was 
calculated as a function of volume fraction (p by means of accelerated Stokesian dynamics simula- 
tions. Numerical results indicate a singular behaviour of the effective viscosity as (p approaches the 
maximal close packing fraction (prep- We quote here from "... The exact form of this singular 

behaviour is not known. Results from lubrication theory for cubic lattices would suggest that the 
singular form should consist both of 1/e and Ine, where e = 1 — {(p/ (prcpY^^ , but the relative amount 
of each term is unknown... As far as we are able to tell at this point, the Ine behaviour accurately 
describes the numerical data." 

When the volume fraction is (approximately) constant over the flow region, e ~ 5. One of the 
objectives of our investigation was to address the issue of the unexpectedly weak blow up in fl^, 
and determine the asymptotic order of the effective viscosity coefficients as 5 — 0. Our analysis of 
the shear viscosity based on the discrete network approximation, showed that ^* = 0{5^'^^'^) 
as (5 ^ 0, while the formal estimate based of the local lubrication analysis would give a higher rate 
0(5~^/^). In dimension three, the corresponding rates, given by the network approximation, are, 
respectively, In 5 and (see |1]). Our analysis offers an explanation of the weak blow up of the 
shear viscosity. This analysis is also consistent with the numerical simulations in Jl] up to the 
difference in the space dimension. 

We also show that the asymptotic order of the extensional viscosity A* depends on the geometry 
of the particle array. For generic disordered arrays the network partitions the domain into poly- 
gons (Delaunay cells), most of which are triangles. We have showed that for these generic arrays 
A* = 0{5-^/^). The same asymptotic rate is obtained for a larger class of networks called quasi- 
triangulated. In a quasi-triangulated network, the percentage of triangular cells may be relatively 
small, but the subnetwork containing triangular cells must be spanning. Another class of networks 
for which A* = 0((5^^/^) consists of rectangular periodic arrays aligned with the boundary of the 
flow. More generally, the same rate is obtained for arrays containing a single spanning chain of 
neighboring particles, perpendicular to the part of the boundary where the velocity is prescribed. 

We show that in case of the strong blow up, the leading term in the asymptotics of A* can be 
uniquely determined by solving a simplified linear system of the network equations, provided the 
array is quasi-triangulated. The simplified system, obtained by neglecting rotational and pairwise 
shear motions of particles, provides an efficient computational tool for evaluating the dependence 
of the effective viscosity on the geometry of the particle array and external boundary conditions. 
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When a network contains neither a spanning triangulated subgraph, nor a spanning vertical path, 
the extensional effective viscosity can be of order This weak blow up is obtained, for instance, 

in the case of a periodic rectangular lattice, with interior edges misaligned with the orientation of 
the prescribed boundary data. 

Our results imply that the ratio of A* to fi* for generic disordered particle arrays is 0{5^^). 
By contrast, in a Newtonian fluid this ratio depends only on the space dimension. Therefore, our 
results indicate possible non-Newtonian behaviour of the effective fluid. 
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A Shear and extensional flows. Ratio of the viscosities in a New- 
tonian fluid 

The following types of flows are relevant to our investigation. 

Shear flow. Consider the steady shear flow a homogeneous fluid characterized by the constant shear 
rate 7. The velocity is given by 

= ( ^0"^ ) . (A.l) 



and the strain rate tensor is 



(A^2) 



Since the stress tensor is symmetric and independent of x, 

E° = ^5°.eO,dx=i5?27|f^|, (A.3) 

where \^\ = jQdx. 

In the case of a homogeneous Newtonian fluid with viscosity fi, = "iiJ-e^i^ — PI, so that 
= ^7- Using the formula (|3.4() we obtain 

= /i, (A.4) 

as expected. 

Extensional flow. In this case the fluid is being extended in the horizontal direction and simulta- 
neously contracted in the vertical direction at the same constant rate e. The velocity is 



ext 



and 



= ( 5 _° J • (A.6) 
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For a homogeneous Newtonian fluid, S^^^ = 2/^6^^^ — PI, and thus S^i = 2/ie — P, S22 = — 2^e — P. 
Next using the definition we obtain 

A* = 4/u. (A.7) 
Therefore, for a Newtonian effective fluid the ratio A*//x* equals 4 (in two dimensions). 

B Proof of Proposition 17.11 

The proposition 17. II wiU be proved if we prove the following 

Proposition B.l Suppose the network graph T is quasi-triangulated. Then every solution of the 
homogeneous system h5.2y\) is of the form twg where t is arbitrary real and wq is given by \5.2^) . 



Proof. First note that (|5.29() is Euler-Lagrange system for the functional 

i ^ 

2 



= • ^ = ^ E E ^%^' - • q'')' + E ■ c^f. (B.l) 



Clearly the minimum of Qhom is zero. Thus every solution of 1)5. 29() is a minimizer of Qhom- On 
the other hand, Qhom{T^^i ■ ■ ■ T^) = if and only if the vectors T* satisfy the system of equations 

(T^ - T^) . q*J- = 0, i = l,2...N, j e K 

• q' = 0, i G /. ^ ' ' 

Therefore, a vector z = (T^, . . . T^)^ solves (HT^ if and only if T% i = 1, . . . , iV solve jH^- The 
solvability of HB.2|) will be directly linked to the geometric structure of the graph F. We begin by 
observing that q' = ±62- Thus the second set of equations in ()B.2|) yields 



T' = fei,iel, (B.3) 

that is, T* are horizontal for all boundary vertices. Next, consider boundary vertices x*,i G /~ 
(these are vertices connected to di}~), and recall that they belong to a path F~, edges of which are 
interior edges of F. Hence, if ii (z I" , then there is at least one ^2 I~ ,12 7^ ii such that x*^ and 
x*2 are connected by an interior edge 6*^*2 _ Using the first set of equations in HB.2|) and HB.3|) we 
obtain 

(T^i - T^2) . q*i^2 = {f^ - f 2)ei • c^''^ = 0. (B.4) 

Furthermore, 6*^*^ is non- vertical, that is, q*i*2 . g-j^ ^ g, which yields i*^ = f^. Since each x*,i € /~ 
is connected to at least one other, we obtain 

T' = tei, ier, (B.5) 

with the same scalar t. 

Since the graph is quasi-triangulated, there exists an interior vertex x'^ £ F , x'^ ^ F~, 
connected to at least two vertices jX^^ ,ii,i2 S /~ by non-collinear edges. Then, using the first 
set of equations in ()B.2|) we obtain 

(T'i-tei).q'i'*i=0, 

(T'l -tei) -q'l'^^ =0. ^ ' 
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Since q'i'*i and q'i'*2 are linearly independent, (|B.6|) yields T'^ = tei. Next, let Gi be the union of 
r~, x'l and all the edges which connect x'^ to r~. Using the definition of the quasi-triangulated 
graph again, we find a vertex x'^, not contained in Gi, and connected to Gi by two non-collinear 
edges. Repeating the argument following (|B.6|) . we see that T'^ = tei. Then we choose G2 to be 
the union of Gi, x'^ and all edges of F which connect them. Repeating the process we find the 
vertex x;^ and continue until we obtain 

T* = tei, i = l,...iV (B.7) 



with the same scalar t. This means that every solution of ()B.2|) is of the form twQ. Since solution 
spaces of ()B.2|) and (|5.29|) are the same, the proposition is proved. 

The following discrete Korn inequality follows immediately from the Proposition IB. 11 

Corollary B.l Suppose that P is quasi-triangulated. Let W C be the one- dimensional sub- 

space spanned by wq from k5.28\) . and let W-^ denote the orthogonal complement of W in 'R?^ . 
Also, let Qhom o.nd A be, respectively the quadratic form defined in \B. 1\) and its matrix. Then 
there is a constant C > such that the Korn-type inequality 

^Az ■ z = Qhomiz) >Gz-z (B.8) 

holds for all z G W-^ . 

Another straightforward corollary is as follows. 

Corollary B.2 Suppose that P is quasi-triangulated. Then the system 

Az = t 

has a unique solution z S W'^ provided i-LW . 
Remark. The projection Pw onto the subspace W is defined by 

„ z • Vi^O 



In terms of vectors T*, 



En rp'i 
^ , . . . , — t=i 1 „, 
r-w[i- , J- ,...± )- — Wo 



Therefore, using the definition of f in terms of R*, we can write the condition f_LTy as 

n N 



J]] = ^ R* • ei = (B.c 



i=l i=l 

The vectors R* in ^^71^ satisfy HB.9|) . so that f in (|^?77|) with defined by ^^71^ is orthogonal 
to W. This gives the unique solvability of the network equations (|f).26|l . 

Corollary B.3 Suppose that P is quasi-triangulated. Then there is a unique z* G W-^ such that 
every solution of h5.26\) is of the form z* + two, where twg G W . 



This means that solution of the network equations 1)5. 26() is unique up to a horizontal translation. 
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C Proof of proposition 17.21 



Proof. First we observe that the form Q is a sum of non-negative terms, each of which corresponds 
to an edge of the network graph T. Removal of an edge from T corresponds to deletion of one 
non-negative term in Q. This means that for each subgraph T' of F, Q(F) > Q(F'). Next we choose 
F' to be the maximal quasi-triangulated subgraph of F. We show that minQ(F') > 0. Indeed, 
min(5(F') = if and only if the corresponding system (|7.1fl)) has a solution. To show that this 
system has no solutions, introduce new unknowns T* = T* — e2. Then from H7.10|) we obtain 

(t^ - T-') . q*^' = 0, (C.l) 



t^.q^ = |-' -h-^^^^ (C.2) 
when i £ / . 



Since the vectors q* are vertical, ()(].2|) yields 

t* = fei, (C.3) 

where f is a constant. Recall that F contains a path F~ which consists of all boundary vertices 
connected to dO,~ and all interior edges connecting these vertices. Hence, each x*i,ii € /~ has 
a neighbor x^^,i2 € I~ and thus (t*^ — t*^)ei • q*i*2 = q from Since two boundary vertices 

cannot be joined by a vertical edge, ei • q*i*2 q_ "Yhis implies that all t^,i G /~ are equal, that is 

f' = tei,i e r, (C.4) 

where i is a constant. Next, consider the boundary path F^. By definition of F' there is a vertex 
connected to the boundary vertices x*i,x*2, ii,i2 G /~ by non-collinear edges of F. 
Then from HC.1|) and HC.4I) we have 

(T-'i -tei) -q^i'*! =0, 
(T-'i -tei) •q^'i'*^ =0. 

Since q-?i'*i and q-?i'*2 are linearly independent, we obtain T-'i = tei. Now this argument can be 
used recursively. Next we choose Gi to be the union of vertices x*,z G I~ ,x^^ and the edges of F' 
which connect these vertices. Repeating the argument, we find a vertex x-^^ ^ q^^ connected to at 
least two vertices of Gi by non-collinear edges, which yields T-^^ = f^^^ a,nd so on, until we obtain 
T* = tei for all vertices x* which belong to F'. By assumption, F' contains at least one vertex 
x^ G I^. But then -q^ = —tei -62 = which contradicts (|(].2|) . This contradiction shows that 
the system (|C.1|) . ()C.2|) has no solutions. 
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